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ABSTRACT 

Direct simulations of wind musical instruments using the compressible Navier Stokes equations have 
recently become possible through the use of parallel computing and through developments in numerical 
methods. As a first demonstration, the flow of air and the generation of musical tones inside a soprano 
recorder are simulated numerically. In addition, physical measurements are made of the acoustic signal 
generated by the recorder at different blowing speeds. The comparison between simulated and physically 
measured behavior is encouraging and points towards ways of improving the simulations. 
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ABSTRACT 

Direct simulations of wind musical instruments using 
the compressible Navier Stokes equations have recently 
become possible through the use of parallel computing 
and through developments in numerical methods. As a 
first demonstration, the flow of air and the generation 
of musical tones inside a soprano recorder are simulated 
numerically. In addition, physical measurements are 
made of the acoustic signal generated by the recorder at 
different blowing speeds. The comparison between sim¬ 
ulated and physically measured behavior is encouraging 
and points towards ways of improving the simulations. 

1. INTRODUCTION 

A flue pipe consists of the flue channel (where the 
jet of air is formed), the labium (sharp edge), and the 
pipe (resonant cavity). The jet of air emerges from 
the flue channel, impinges the labium, and oscil¬ 
lates around the labium generating acoustic waves 
which travel into the pipe and return to interact with 
the jet of air in a complex feedback cycle. Predict¬ 
ing the acoustic frequencies which are generated 
by a flue pipe is generally difficult, and is possible 
only in simple cases and well-studied geometries. 
The study of flue pipes dates back to Helmholtz 
and Rayleigh and perhaps earlier, and remains very 
active today because the problem is very difficult. 
Approximate physical models of flue pipes have 
been developed over the years which work well 
in certain cases, but not always, and the assump¬ 
tions behind the models are not fully understood 
(Hirschberg [1], Verge [2] and references therein). 

In this paper, a soprano recorder is simulated using 
the compressible Navier Stokes equations. As far as 
we know, these are the first simulations of flue pipes 
using the compressible Navier Stokes equations to 


be presented (also see Hirschberg [1]). Further, the 
predictions of the simulations are compared against 
physical measurements of the acoustic signal gen¬ 
erated by the recorder at different blowing speeds. 
The comparison is encouraging and points towards 
ways of improving the simulations. 

Some of the difficulties of simulating subsonic hy¬ 
drodynamics with acoustic waves are summarized 
in Hardin [3]. Most importantly, one must choose 
the time step of integration very small in order to 
follow the fast-moving acoustic waves. Second, 
one must control the slow-growing numerical in¬ 
stabilities which arise in simulations of subsonic 
flow. Third, one must impose boundary condi¬ 
tions around the simulation region that do not reflect 
acoustic waves. 

The requirement of small time steps in simulating 
subsonic flow makes the computation expensive, 
but it also simplifies the computation. In particular, 
explicit numerical methods can be used which are 
simple to program and ideal for parallel comput¬ 
ing. In many other PDE-based problems, explicit 
methods are usually avoided because they require 
very small time steps to remain numerically stable. 
In the present problem however, the time step must 
be kept very small anyway to follow the acoustic 
waves, so there is a match between the problem and 
the use of explicit methods. 

The remaining parts of this paper are as follows. 
First, we describe the physical model and the ge¬ 
ometry of the soprano recorder which we have sim¬ 
ulated numerically. Then, we outline our numerical 
methods. Subsequently, we present representative 
results from the simulations of the recorder, and 
a comparison with physical measurements of the 
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acoustic signal generated by the recorder. Finally, in the discrete numerical system (see [4]). They are 
we comment on our results and on future directions. used to prevent the reflection of acoustic waves. 


2. PHYSICAL MODEL 

The flow of air is modeled using the adiabatic ap¬ 
proximation [4] so that the fluctuations of pressure 
are proportional to the fluctuations of density times 
the square of the speed of sound. The density p 
and the components of the flow velocity in the x,y 
directions 14, K, describe the flow (/>. \. .. I . are 
functions of space and time). The following Navier 
Stokes equations are used, 


Regarding the initial blowing of air into the flue 
channel, it is chosen to rise smoothly to a final 
steady value. The following formula is used both 
for the velocity and the density at the inlet and the 
outlet. 
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The constants v = 0.15 cm 2 /s and c s = 34400 cm/s 
are the kinematic viscosity and the speed of sound 
corresponding to air flow at 22 degrees centigrade. 


3. NUMERICAL METHODS 

The Navier Stokes equations are solved numerically 
using a uniform grid and two different explicit nu¬ 
merical methods: the lattice Boltzmann method and 
a traditional finite difference method [4, 5]. The 
simulation results obtained with the lattice Boltz¬ 
mann method and the finite difference method are 
very similar to each other; and for the sake of 
brevity, only the lattice Boltzmann simulations are 
presented here. 


The simulations are based on the two-dimensional 
model of the soprano recorder which is shown in 
figure 7. A steady supply of air is forced through 
the flue channel and exits through the outlet at the 
top of figure 7. The walls around the flue pipe (gray 
areas in figure 7) are non-slip walls which reflect 
acoustic waves. By contrast, the outlet and the two 
vertical “walls” above the flue pipe do not reflect 
acoustic waves so that the region above the flue pipe 
resembles an infinite space as much as possible and 
not another resonant cavity. 

Non-reflecting boundary conditions are imple¬ 
mented by imposing both the density and the ve¬ 
locity at the boundary. For example, ambient pres¬ 
sure is imposed at the outlet, and a pressure which 
is slightly larger than the pressure drop of Hagen- 
Poiseuille flow through the narrow flue channel, 
is imposed at the inlet. Further, a parabolic ve¬ 
locity profile is imposed both at the outlet and 
the inlet. These boundary conditions are overcon¬ 
strained; however, they do not cause any problems 


A fourth-order artificial viscosity filter is used in 
combination with the numerical methods in order to 
control slow-growing nonlinear instabilities [4, 6]. 
The filter does not alter the fluid dynamics at long 
wavelengths, but provides a small amount of ad¬ 
ditional dissipation at short wavelengths compara¬ 
ble to the grid spacing. Similar filters have been 
very common in transonic and supersonic flow 
(Peyret&Taylor [7]), and appear to be necessary 
for simulating compressible flow in general. 

In the present simulations, the spatial resolution 
is Ax = 0.01 cm so that 10 fluid nodes fit along 
the width of the flue channel. About 0.8 million 
fluid nodes are used in the simulation, and they are 
computed in parallel by 22 non-dedicated Hewlett- 
Packard HP9000/715 workstations (see [6, 8]). The 
time step is At = 2.17 x 10 -7 s so that about 150000 
integration steps are needed to produce about 32 ms 
of simulated time. Each integration step takes about 
1 s of runtime, and a typical simulation of 150000 
integration steps lasts approximately 2 days. 
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4. RESULTS 

Figure 7 shows the 2D geometry of the soprano 
recorder which is used in the simulations. The real 
recorder is cylindrical and looks like figure 7 when 
sliced vertically [4]. The recorder is 20 cm long and 
is closed at the far end. Thus, the recorder behaves 
as an open-closed pipe with major frequencies in 
odd ratios 1:3:5 etc. The two lowest modes of an 
ideal open-closed pipe 20 cm long are 430Hz and 
1290Hz, and the present recorder produces frequen¬ 
cies around 400Hz and 1100Hz as can be seen from 
the simulations and the physical measurements be¬ 
low. 

Figures 1 to 6 show acoustic signals obtained from 
simulations and from physical measurements at dif¬ 
ferent blowing speeds. The figures are matched in 
pairs (1&2, 3&4, 5&6) based on the resemblance 
between the simulated and the physical signals. A 
matching between simulated and physical signals 
based on the blowing speed is not possible prob¬ 
ably because the simulations are two-dimensional, 
and also because the blowing speed in the physical 
setup is only known very roughly. Pair 1&2 shows 
the behavior at low blowing speed where the funda¬ 
mental mode 400Hz dominates. Pair 3&4 shows the 
behavior at medium blowing speed where the fun¬ 
damental 400Hz and the “third” 1100Hz have ap¬ 
proximately equal size. Finally, pair 5&6 shows the 
behavior at high blowing speed where the “third” 
1100Hz dominates. 

5. CONCLUSION 

We have presented representative results from the 
simulation of a soprano recorder using the com¬ 
pressible Navier Stokes equations. Physical mea¬ 
surements of the acoustic signal generated by the 
recorder at different blowing speeds are in reason¬ 
able agreement with the predictions of the simu¬ 
lations. However, there are also important dif¬ 
ferences between the simulated and the physical 
behavior. One source of differences is that the sim¬ 
ulations are the initial transient behavior, and the 
physical measurements are the steady-state behav¬ 
ior. The present simulations can not be continued 
much longer than shown here because vorticity ac¬ 


cumulates above the flue and affects the operation 
of the flue pipe. In the physical world, the gen¬ 
erated vorticity quickly moves away from the flue 
because of air entrainment and also because of the 
available volume along the third dimension. 

In figure 3 (medium blowing speed) we can see that 
good acoustic oscillations are only achieved for a 
few cycles, and then are lost. We suspect that this is 
a result of the accumulation of vorticity above the 
flue because the presence of vortices changes the 
ambient pressure level outside the flue pipe. Work 
is in progress to find ways of removing the gener¬ 
ated vorticity, and to understand better the simula¬ 
tion results with regard to boundary conditions and 
the two-dimensionality of the simulations versus 
the three-dimensionality of the physical world. 
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Figure 1: Simulation, 811 cm/s 
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Figure 3: Simulation, 953 cm/s 
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Figure 2: Physical, 734 cm/s 
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Figure 4: Physical, 1140 cm/s 
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Figure 5: Simulation, 1108 cm/s 
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Figure 6: Physical, 1558 cm/s 


The above figures (1 to 6) show acoustic signals obtained from simulations and from physical measurements 
of the 20 cm closed-end recorder at different blowing speeds (mean flow speed through the flue channel). 
The simulation data is sampled 4 cm above the labium, and the physical data is sampled 100 cm above 
the labium. Each figure plots the time series of the normalized density fluctuation at the bottom, and the 
frequency spectrum (decibels with respect to standard pressure level) at the top. The physical measurements 
have arbitrary units of amplitude because the measurements were not calibrated. 
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Figure 7: The two-dimensional geometry of the soprano recorder which is used in the simulations. The 
numbers shown correspond to millimeters. The pipe is located at the bottom of the picture, and measures 
20 cm long and 1.34 cm wide. The flue channel is located at the bottom left corner, and measures 4 cm 
long and 0.1 cm wide. At a distance of 0.4 cm in front of the flue channel, there is a sharp labium which 
measures an angle of 14 degrees approximately, and is positioned slightly below the midline of the flue 
channel. Namely, the tip of the labium is located at 1.34 cm from the bottom of the pipe, and the flue 
channel is located between 1.3 cm and 1.4 cm. 



Figure 8: Simulation of a 20 cm closed-end soprano recorder. Picture is 30 ms after startup. Iso-vorticity 
contours are plotted. The dashed lines show the decomposition of the problem into 22 workstations. The 
gray-shaded areas are not simulated. 
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